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Nonlinear evolution of one-dimensional planar perturbations in an optically thin radiatively 
cooling medium in the long-wavelength limit is studied numerically. The accepted cooling function 
generates in thermal equilibrium a bistable equation of state P(p). The unperturbed state is taken 
close to the upper (low-density) unstable state with infinite compressibility (dP/dp — 0). The 
evolution is shown to proceed in three different stages. At first stage, pressure and density set in 
the equilibrium equation of state, and velocity profile steepens gradually as in case of pressure-free 
flows. At second stage, those regions of the flow where anomalous pressure (i.e., with negative 
compressibility) holds, create velocity profile more sharp than in pressure-free case, which in turn 
results in formation of a very narrow (short-wavelength) region where gas separates the equilibrium 
equation of state and pressure equilibrium sets in rapidly. On this stage, variation in pressure 
between narrow dense region and extended environment does not exceed more than 0.01 of the 
unperturbed value. On third stage, gas in the short-wavelength region reaches the second (high- 
density) stable state, and pressure balance establishes through the flow with pressure equal to 
the one in the unperturbed state. In external (long-wavelength) regions, gas forms slow isobaric 
inflow toward the short-wavelength layer. The duration of these stages decreases when the ratio 
of the acoustic time to the radiative cooling time increases. Limits in which nonlinear evolution of 
thermally unstable long-wavelength perturbations develops in isobaric regime are obtained. 
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I. INTRODUCTION 



Thermal instability is thought to be a mechanism which causes a radiatively cooling medium (RCM) to break onto 
two phases: relatively dense cold filaments and clumps surrounded by rare hot plasma. This phenomenon is common 
for both cosmic and laboratory optically thin plasmas which lose their energy radiatively (see for more discussion 
A comprehensive linear theory of thermal instability has been developed in the pioneering paper by Field The 
instability results from the fact that the radiative energy loss rate grows with density. A local increase of the density 
brought in a uniform equilibrium plasma, enhances cooling rate and stimulate temperature to fall. This leads in turn 
to subsequent progressive growth of the density in order to keep pressure constant. 

Last decade, nonlinear aspects of thermal instability have been extensively investigated both analytically and 
numerically, and considerable progress in understanding of the dynamics of phase transitions and structure formation 
in bistable media was achieved. 

The equations governing dynamics of RCM, besides the nonlinearities intrinsically connected with gas dynamics 
processes, contain also nonlinearities related to radiative energy loss rate, E(p, T), which is square function of density 
and, as a rule, strongly nonlinear function of temperature. In general, nonlinear dynamics of such media is complex 
for analytical study, and several approximations were used to reduce the equations to more simple and useful form. 
This procedure is based on the presence of different characteristic time scales in RCM: characteristic radiative cooling 
time tr, conductive time r x , and acoustic time ta- When these times differ considerably from each other one can find 
the fastest process which relaxes to its equilibrium state, and thus eliminate one of the governing equations. In this 
manner Doroshkevich and Zeldovich ||, assuming the medium to follow isobaric dynamical regime, have described 
dynamics of RCM by a nonlinear heat-diffusion equation 



dT d dT 

which belongs to a class of diffusion-reaction equations, here the reduced cooling function e(T) represents the "reac- 
tion" term. In subsequent papers Meerson and Sasorov Q , and then Meerson et al. j5| have developed this approach to 
reduce the dynamics of RCM to more simple set of equations depending on interrelation between characteristic times. 
In the short-wavelength limit when acoustic time is the shortest, the dynamical regime is isobaric and the equations of 
motions are reduced to a diffusion-reaction equation in Lagrangian variables |^| , and contrary, when acoustic time is 
the longest (the long- wavelength limit), the dynamics of RCM is described by equations of gas dynamics with equation 
of state determined by "thermal equilibrium", i.e., by balance between radiative cooling and external heating [Q. 

The above approximations suggest that the wavelength of perturbations of hydrodynamical variables satisfies one 
of the two strong inequalities: ta *C tr for the short-wavelength and ta 3> tr for the long-wavelength limits, 
respectively. Perturbations with wavelengths close to the characteristic length, i.e., ta ^, tr or T4 tr, show complex 
behavior and require in general numerical study treated with the full set of equations. We address in this paper 
the question of what are the wavelength limits in which isobarical approximation can be used for description of 
thermal instability. We consider perturbations which initially are small (i.e., linear), and show that these limits 
are rather broad: at t > 3 — 15 ta, perturbations with wavelengths much larger than wavelengths determined by 
Meerson et al. [|J, i-e., ta ^S> tr, evolve isobarically. Our conclusion is therefore that at comparatively late stages 
long-wavelength perturbations can be described by a diffusion-reaction type equation for temperature and density, 
equivalent to the equation for isobaric short- wavelength perturbations. Initial stages of long- wavelength perturbations 
can be understood in terms of wave dynamics in a pressure-free medium, and are reasonably described by nonlinear 
advectivc-diffusion equation for velocity 

The paper is organized as follows: in Section II we present the basic equations and discuss limiting cases which 
allow analytical description, in Section III we show numerical results for long-wavelength perturbations. Section IV 
contains a summary. 



II. GOVERNING EQUATIONS, CHARACTERISTIC TIMES, LIMITS 



A. Equations 

The basic equations governing dynamics of a radiatively cooling gas were proposed originally by Field and include 
in the energy equation source and sink terms describing heating by external sources and radiative energy losses. In 
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planar geometry the equations are written 



(1) 



.' du du\ dP d 2 u 
1 fdP dP\ j du , d ( dT\ ( du\ 2 

P - —pT = 0, (4) 
M 

where C is the net cooling function defined as C(p, T) — pA(T) — T(p, T), with T, as a rule, weekly dependent on p 
and T, other notations have a common meaning. 



B. Characteristic times 



Three characteristic time scales are relevant for description of perturbations in RCM: acoustic time ta = A/cq, 
conductive time t x = IZpX 2 / p\, and radiative cooling time tr = IZT/App, here cq = \J jPq / po is the adiabatic sound 
speed. It is obvious that 

t x ~ 3A 2 /£c , 

where I is a free-path of gas particles, and thus in framework of hydrodynamical description, i.e., A>l, 



T x = \T A /£ > T A . 

For cooling rate per unit mass pA(T) ~ (TiCoAEn/ p, where Oi is the cross-section of inelastic collisions, AE, the 
energy lost by particles in a single inelastic collision, n is the gas number density, characteristic cooling time is 



TIT i 
eAE Co 

where e = Ui/a e < 1, the ratio of inelastic-to-elastic cross-sections. Since energy transfered in a single collision, as a 
rule, is small compared to thermal energy AE < 1ZT, we get from here ^ £/cq. 



C. Limits 



Therefore, three different intervals of wavelengths are relevant, see Fig. 1: short wavelength s A# < A < A^, 
intermediate wavelengths Af < A < Xa, and long wavelengths A > A^, here Af = \ZxoTq/ pqAq, the Field length 
determined as a length where conductive heat transfer and radiative losses equate, i.e., tr = r x , A^i = t^cq is the 
acoustic length where ta = tr and Xh ~ £ is the wavelength at which hydrodynamical description violates. It is 
seen from Fig. 1 that always the interrelation between characteristic lengths Xh < Xp < Xa is valid. From this point 
of view, the "direct crossover" regime corresponding to the case of large heat conduction, A^ <C Xp, is out of the 
hydrodynamical approximation. 

The principal difference between the short and intermediate- wavelength intervals is in more efficient conductive heat 
transfer compared to radiative losses for short- wavelength perturbations, which leads to erasing of a monochromatic 
perturbation (a normal mode with fixed wavenumber k = 2ir/X, e.g., ST oc e lkx ) with A < Af (see [0j2|). A non- 
monochromatic localized isobaric perturbation (a wave packet of a set of normal modes, e.g., ST oc J STke lkx dk) with 
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characteristic size L < Xp shows decrease in amplitude at initial time, which then changes to increase due to thermal 
instability of long- wavelength modes presenting in the perturbation Q. Since in these intervals ta is the shortest 
characteristic time, pressure balance sets in first over the perturbation equations reduce to the isobaric form [Bp) 



dT 8 f-dT\ 

dt ^ dm \ " dm J 



(5) 



where 



L{T,P) 



7~ 1 M 
7 K 



is the reduced cooling function, 



- = 7 - 1 X 
X 7 PT 1 



the reduced heat conduction coefficient, Tdm — dx, P =const. 

In the long- wavelength limit, acoustic time ta is much larger than characteristic radiative cooling time tr, and 
conductive heat transfer is negligible: t x ~ ) 2 /x » T ,4 T R, see Fig. 1. In this limit the equations can be reduced to 
more simple form describing gaseous medium with a given (anomalous) equation of state M . Indeed, when ta ^> tr 
the energy balance sets in on the shortest time tr, and thus temperature and pressure relax rapidly to the values 
determined from the equation of steady energy balance 



Note, that Eq. (6) with cooling function A(T) shown in Fig. 2 and T =const maps (A, T) plane on (P, p) plane 
to give the equilibrium equation of state P — P{p) of a van der Waals type. Equilibria with temperature in the 
interval where dhi K{T)/d\nT < 1 are thermally unstable, and correspond to anomalous part of the P(p)-curve with 
dPj dp < [0j2||] . These arguments allowed Meerson and Sasorov j| to reduce the momentum equation of RCM in 
the long- wavelength limit to the one with a given ("equilibrium") equation of state P = P{p) 



where P = P(p) is determined from Eq. (6). For cooling rate with din A(T)/dlnT > 1 below the temperature T = Tl 
and above T = Tu, and <iln A(T)/dlnT < 1 in the intermediate interval, the equation of state has a form similar 
to the van der Waals equation with two rising (stable) branches of the curve P{p) in the low density p < pu (high 
temperature T > Tjj) and high density p > pl (low temperature T < Tl) ranges, and a falling (unstable) branch of 
the P(p)-curve in the intermediate range pu < P < Pl {Tl < T < Tu). A^-shaped curves in Figures 4 and 6 show the 
"equilibrium" equation of state P(p) corresponding to Eq. (6) with the cooling function A(T) depicted in Fig. 2 (see, 
also [0j6|). This equation of state is commonly used for description of pattern formation in a bistable medium in the 
long- wavelength limit [Q. (For the interstellar medium T L - 10 2 K and T v - 10 4 K, see HfuJ). 

We will mostly concentrate in this paper on the dynamics of perturbations in an RCM which is in the upper 
marginally stable point: T ~ Tu, dlnA/dlnT ~ 1, thus even small (linear) perturbations break the initial state. In 
this state an RCM can be prepared, e.g., by adiabatically slow decrease of the intensity of external heating source 
in an RCM kept at fixed pressure Pu, or adiabatically slow increase of pressure at fixed intensity of heating sources. 
As at the equilibrium point T = Tu the compressibility is infinite, i.e., dP/dp = 0, small perturbations around this 
point evolve to form localized structures as was demonstrated by Meerson and Sasorov Q . In the next order positive 
perturbations of density decrease pressure, providing anomalous (negative) compressibility 



where 5p — p — pu, which amplifies the nonlinearity connected with the advective term in momentum equation (3). 
It is clear that nonlinear term in Eq. (8) does not saturate the instability, and therefore both the advective term and 



C(p,T) = 0. 



(6) 




(7) 



P( P ) = Pu- Y 6p 2 , 



(8) 
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anomalous negative compressibility can be stabilized only by a strong nonlinearity in the equation of state P(p), which 
gives rise to the second stable branch. This saturated state corresponds therefore to a new state of a medium separated 
onto two phases. The difference between densities of two phases being in pressure equilibrium is considerable, as a 
rule about two orders of magnitude, and only numerical study can be applied to these strongly nonlinear saturation 
stages H . In order to study dynamics of long- wavelength localised structures in a fluid with negative compressibility 
Meerson and Sasorov j| have approximated the falling branch of the equation of state as P(p) = Pjj + const/p, and 
described the so-called unlimite d instab ility 

Initial stages, when \Sp\ <C \J1P\j /a 2 , can be described by the advective-diffusion equation 



( du du\ d 2 u 

Numerical simulations show that equation (9) gives qualitatively reasonable description of velocity profiles up to 
t £j (2 — 15)ta depending on wavelength. Within this time interval pressure varies from Pjj to minimal value P = 
Pu — AP m , and arised pressure gradients force then velocity to steepen sharply. It forms, in turn, a short-wavelength 
region, where thermal regime separates the "equilibrium" equation of state P = P(p) and isobaric distribution sets 
in rapidly, and then velocity relaxes gradually to profiles adequate to isobaric dynamics. We argue therefore, that at 
larger times long-wavelength perturbations in RCM evolve isobarically, and can be described by equation (5). 



III. NUMERICAL RESULTS 



To confirm this conclusion we show here two numerical models which follow the evolution of perturbations from 
initial linear state to final asymptotic behaviour. In numerical simulations we have used a 6-parametric finite difference 
scheme of Eqs. (l)-(4) which is free of nonphysical effects of numerical viscosity and absolutely stable (the detailed 
description of the scheme will be given elsewhere). 

To make presentation more clear we use dimensionless variables with length, time and velocity normalized, respec- 
tively, to the initial size A of a perturbation (which is assumed to be localized), to the dynamical time ta — A/co, and 
to the adiabatic sound speed cq; hydrodynamical variables p, P and T are normalized to their unperturbed values. 
In these units cooling time is equal to fa — tr/ta- As we concern here dynamics of long- wavelength perturbations, 
the characteristic size is A > Xp, (in numerical examples given in the paper Xp = 0.001A). Dimensionless cooling 
function A(T) has been taken in the form 



A(T) = {ai + a 2 tan[|(T*-T)]} , (10) 

as shown in Fig. 2. This function represents typical cooling function for bistable media, and generates equation of 
state of the van der Waals type when "thermal equilibrium", i.e., pA(T) = T(p,T) is reached; we assume through 
the paper T(ft, T) = A(l)=const. The unperturbed state was taken near the upper unstable point. We present here 
results with initial conditions taken in the form of small velocity perturbations corresponding to a localized symmetric 
inflow to the center of a grid zone 



, s _ J -n sin[fe(a; - x )], \x - x | < 27r/fc, , , 

u[X,U)-<^ ^ otherwise, [ > 

while perturbations of p, T, and P were zero. This choice of initial conditions seems to be most appropriate in 
study of long- wavelength perturbations. Due to radiative losses initial conditions corresponding to perturbations of 
density, temperature and pressure are inevitably more affected by the procedure of preparation such perturbations. 
We therefore, do not consider here alternative dynamics which would start from small perturbations of density, 
temperature and pressure. 

For adopted type of initial conditions velocity perturbations disturb density, temperature and pressure, and radiative 
cooling drives them to relax in a short time to values determined by thermal equilibrium (6). Arised gradients of 
pressure redistribute then velocities, such that the resultant velocity field tends to set the flow in isobaric regime. At 
late stages the flow is approximately isobaric through the perturbed region. The boundary conditions used correspond 
to a "closed" box, i.e., zero velocity and gradients of density and pressure. In most models the grid zone size was 
20 times of the size of a localized perturbation, however test runs with larger size (80 times of the perturbation) 



5 



show results to be insensitive to the grid size. This is a direct consequence of the fact that the compressibility in the 
unperturbed state (which is taken to be close to the upper unstable point) is very high, and a localized perturbation 
propagates slowly: the acoustic velocity cu = y/ (dP/dp)u ~ 0.25. Moreover, the initial perturbation brings most 
fraction of perturbed gas into a state with negative compressibility where the perturbation does not propagate. The 
time required for perturbations to reach a boundary is therefore at least 80 dynamical times, which is larger than 
the total computational time in all cases we considered. This means, that the boundary conditions we used weakly 
affect dynamics of perturbations. (Note, that systems whose unperturbed states lie on the low-density stable branch 
well below the upper unstable point, can be brought into thermal instability only by perturbations of finite, i.e., not 
small, amplitude. We do not consider such cases in this paper.) 

Fig. 3 shows dynamics of a perturbation with initial amplitude of velocity uq = 0.03 and with tr = 0.33 (i.e., A = 
3\a), thus it corresponds to the long-wavelength limit. On the upper panel velocity profiles are plotted. At initial 
times, t < 10, as gas in central regions is involved in motion the velocity amplitude decreases. In gases with normal 
equation of state, i.e., dP/dp > 0, perturbation decays onto two acoustic waves propagating in opposite directions, 
and a standing and damping entropy mode [in case of perturbations given by eq. (11), decayed waves represent 
two pieces of sinusoidal waves with velocity amplitude of uq/2]. Radiative losses with tj; « 1 modify the inflow, 
transfcring the internal (compressed) gas onto the anomalous part of equation of state (dP/dp < 0), and thus make 
the perturbation locally growing, i.e., non-decaying. At times 10 < t < 15 the velocity profile steepens and motion 
concentrates in a rather narrow region with increasing density, as shown in the bottom panel. 

Relaxation of pressure to the equilibrium (anomalous) equation of state P(p) occurs at times t ~ tr (these stages 
are not shown). Increasing density in central regions amplifies pressure gradient through the perturbation, and thus 
forces the velocity amplitude to increase. It results in more sharp steepening of velocity profiles compared to the 
pressure-free case as seen in Fig. 3a for time interval 10 < t < 15: the profile for t — 15 clearly shows two inflection 
points outside the central region very close to peak values, which is not what can be expected for pressure-free motions, 
or motions with normal equation of state. In turn, the inflow with increasing velocity amplitude forms very narrow 
region (a "droplet" with sizes Ax <C A, as seen in Fig. 3b), through which pressure equilibrium sets in rapidly. As 
a result, at these stages structure of the perturbation can be understood as distinguished onto two regions: first, 
the "long-wavelength" part outside the "droplet" where the equilibrium equation of state dominates, and second, the 
"short-wavelength" part (the "droplet" itself) where motion tends to set in the pressure equilibrium. The two regimes 
manifest themselves clearly in Fig. 4, where loci P — p connecting pressure and density through the computational 
zone for different times, i.e., P(x,t) and p(x,t) profiles with eliminated x, are plotted: the short piece of line at 
the upper equilibrium point shows P — p locus for t — 10. The locus corresponding to t = 15 is readily seen to 
consist of two parts: a small piece of a P — p curve starting from the unperturbed state with decreasing pressure and 
density almost parallel to the growing branch of the equation of state (corresponding to external, long-wavelength 
regions of the perturbation where density and pressure decrease), which then turns around and goes with increasing 
density and almost constant pressure (this part corresponds to internal, short- wavelength regions). Thus, at these 
stages pressure separates the equilibrium curve P(p), and subsequent evolution is characterized by settling the overall 
flow down to isobaric regime. At t Si 20 most fraction of density perturbation concentrates in the isobaric "droplet" 
with characteristic size of \x — xq\ = 0.2 and in the conductive interface, which are in pressure equilibrium with gas 
outside (variation in pressure through the perturbed region is less than 0.001, and main contribution comes from 
surrounding regions where density and pressure lie on the equilibrium equation of state). (Strongly speaking, in 
models with "closed" box boundary conditions we used, pressure depends on time. However, as mentioned above, in 
all considered cases asymptotic behaviour of perturbations establishes on time interval shorter than the time required 
for a perturbation to reach boundaries. One can think, therefore, that isobaric regimes asymptotically established in 
numerical models with "closed" box boundaries give qualitatively correct description of the situation in systems with 
free boundary and fixed external pressure.) 

After the short-wavelength droplet is formed, velocity decreases to u < 0.01 with profiles corresponding to a weak 
inflow, so more gas continues to increase mass of the droplet. Subsequent evolution follows the isobaric regime with 
P — p loci practically coincident with that reached at t = 20, as shown in Fig. 4. Note that in models with fixed 
equilibrium equation of state P = P(p), see jl],^), the established inflow velocity is, in general, higher than in our 
calculations with energy equation (3) solved explicitly. This is due to the fact that in models with the equilibrium 
equation of state, P = P{p), pressure reaches minimal value in the interface layer, and thus at equal conditions 
(e.g., viscosity and conductivity) arised pressure gradients generate higher acceleration. 

In Fig. 5 and 6 we show results for a perturbation with larger spatial size: A = 30A.4, i.e., tr = 0.033. Initial 
conditions are same as in the previous case. The perturbation is seen to reveal similar behavior as that of smaller 
size, however, the dynamics is more violent. The difference is mainly in time scales and the amplitude of velocity: 
perturbations of longer wavelength develop higher velocities and pass same stages quicker. The enhanced steepening 
of the velocity profile due to negative compressibility occurs at t £ 2. Here the difference between enhanced and 
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normal steepening is more obvious than in previous case: rather flat increase in velocity in outer regions, and very 
sharp growth near the symmetry plane. At t ~ 3 the inflow with amplitude u ~ 0.08 forms a narrow short-wavelength 
droplet with almost constant pressure inside [as seen from P — p loci on (P, p) plane, Fig. 6]. However, as the velocity 
amplitude developed is large, pressure in the droplet is higher than the unperturbed value by AP ~ pu 2 ~ 0.01, 
and in very central parts of the droplet overpressure reaches the value AP ~ 0.03 due to higher densities [note, that 
P — p locus in the central overpressured region goes along the equilibrium equation of state P(p) corresponding to the 
second stable brach]. On next stages, 3 < t £ 4, the overpressure leads to expansion of the short-wavelength region 
(with almost isobaric pressure distribution inside) to relax to the equilibrium with external regions. This relaxation 
occurs more violent than in the previous case, due to larger variations in pressure, AP| ~ 0.01 — 0.03 around the 
unperturbed state. As a result, subsequent growth of gas mass in the short-wavelength (droplet) region proceeds in 
a damped oscillatory regime, as seen on (P, p) plane in Fig. 6. 

The dependence of dynamics of perturbations on wavelength (or equivalently on tr) can be understood qualitatively 
from simple estimations of steepening of velocity profiles. The characteristic formation time of a short-wavelength 
region (a droplet) due to steepening can be estimated from at 2 , ~ 2A, where a is the acceleration caused by the 
negative compressibility, 



rvpj AP m 

a- ~ — , (12) 

Po Po* 

AP m is the maximal deviation of pressure from the unperturbed value due to density increase. This gives in dimen- 
sionlcss units 



The deviation of pressure from the unperturbed value increases along the equilibrium P(p)-curve until radiative cooling 
time tr is less than the acoustic time for the shortest wavelength A s generated due to steepening. Thus, AP m can be 
determined from the condition tr = A s . To estimate X s we assume that perturbation of density caused by initial inflow 
(11) in dimcnsionless units is 6po(x) ~ Uq(x). Then, due to thermal instability it grows as 6p(x,t) ~ uq exp(t/r^j). 
On the other hand, one can expect that \ s ~ Aexp(— t/rR), thus giving Sp ~ u /X s . Substituting here AP m ~ c 2 A Sp, 
where c 2 A — \dP/dp\ estimated on the anomalous branch of equation of state [c 2 A <~ 0.1 at p ~ 1.5] 



AP m = ^, (14) 



We obtain therefore for td- 



U - J^r, (15) 

which provides with reasonable accuracy an estimate for the time when the numerical solutions turn from the regime 
with the equilibrium equation of state, P(p), to the isobaric regime. 

It is clear from these estimates and from numerical results presented, that at small values of the ratio uqc a /tr 
(less than 0.1) long-wavelength pertubations develop quiescently as isobaric. In dimensional units the condition for 
long-wavelength perturbations to evolve in a quiescent isobaric regime is written as 



Z±<o.iBl(Bl) . (16) 

TR u \c A J 

When UqC 2 a /tr > 0.1, high pressure gradients developed at initial "equilibrium" stages [i.e., along P(p) equation 
of state] give rise to very fast steepening of velocity profiles, and as a result to formation of the overdense and 
overpressured central short- wavelength droplet, which then relaxes to the equilibrium in an oscillatory regime with 
variations in pressure of |AP| ^ 0.03 around the unperturbed value. Regimes with the equilibrium equation of state, 
P(p), can exist for long- wavelength perturbations only at initial stages, t ^ td- Subsequent evolution implies generation 
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of a short- wavelength core, which in turn results in separation of pressure of the equilibrium curve P(p). Formation 
of such a short-wavelength region is an inevitable consequence of the "unlimited" instability when narrow localized 
structure develops in a medium on the anomalous branch of the equilibrium equation of state 4 . In Fig. 7 we show 
the dependence of pressure in the central droplet on time for several values of tr and fixed c\ — 0.1 at p = 1.5. 
The amplitude of pressure oscillations is clearly seen to increase for smaller tr: approximately, at initial times, when 
the central droplet follows the equilibrium equation of state, AP ~ AP m <~ u c 2 a /tr. The period of oscillations is 

T ~ X s ^pd/AP, where pd — 5 is the density on the second stable branch; at initial times T ~ 2.24f^/ 2 / ^u c\, 
however then it grows due to damping of AP. 



IV. SUMMARY 



We have considered nonlinear evolution of initially small (linear) perturbations in a radiatively cooling medium in 
the long-wavelength limit, i.e., when acoustic time is larger than radiative cooling time. We have not been assuming 
pressure and density to satisfy the equilibrium equation of state, but studied the problem numerically with the 
radiative energy loss term in energy equation treated explicitly. For initial conditions in the form of small localized 
velocity perturbation around a uniform equilibrium close to the upper unstable point, we have shown that early 
stages of evolution tend to proceed along the equlibrium equation of state, P{p), with velocity profiles to steepen as in 
pressure-free case. At intermediate stages density increases enough and anomalous (negative) compressibility becomes 
important and results in more sharp steepening of velocity profiles, which in turn forms a narrow (short-wavelength) 
region where pressure separates the value determined by the equilibrium equation of state and sets in approximate 
isobaric distribution rapidly. One may say that thermal instability works as a regulator which prevents large pressure 
gradient to arise. The maximal amplitude of pressure developed in the perturbation when it separates the equilibrium 
curve P(p), is proportional to the amplitude of initial velocity perturbation and inversely proportional to the radiative 
cooling time [as described by Eq. (14)] 



AP m ~ P ^ T -± (2A) . (17) 

CO TR \ C J 

At late stages, gas in the short-wavelength region reaches the high-density (stable) equilibrium branch of the equation 
of state, and subsequent evolution corresponds to the isobaric inflow of low-density gas onto the high-density core 
(the droplet). We have shown thus, that nonlinear evolution of long- wavelength perturbations after the time interval 



td ~ ~~~ \l 2tatr—— , (18) 

CA V "0 



reaches quiescent isobaric regime, if the wavelength is restricted from above as 



A <0.1A^f^V. (19) 
uo \caJ 

In these limits, long-wavelength perturbations in thermally unstable media can be described as isobaric by the reduced 
equation (5). Perturbations with wavelengths longer than the limit (19), after t — td when pressure separates the 
equilibrium equation of state, develop localized oscillating structures with initial pressure amplitude determined by 
Eq. (14) and damped then due to viscosity. 
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FIG. 1. The dependence of characteristic times on wavelength. Point F corresponds to the Field length X F , A — to the acous- 
tic thermal length Aa, H — to the length where hydrodinamical description violates \h ~ (■■ The interval Ajj < A < \f corre- 
sponds to the short-wavelength limit, Af < A < Aa — to the intermediate wavelengths, and Aa < A — to the long-wavelengths. 

FIG. 2. The cooling function A(T) used in the paper: ai = 0.915, Q2 = 0.61, T* = 1.11. The point indicates the unperturbed 
state; the heating rate was taken to balance energy losses in this point. 

FIG. 3. Nonlinear evolution of perturbations in RCM: A = 3Aa, tr = 0.33. Upper panel (a) shows velocity profiles for the 
times t = 0, 10, 15, 20. The amplitude of initial velocity perturbation is uo = 0.03. Density, temperature and pressure have 
not been perturbed. Bottom panel (b) shows density profiles corresponding to the velocity profiles given on the upper panel 
[initial density profile po(aO=const is not shown]. 

FIG. 4. P — p loci [obtained from P(x,t) and p(x,t) profiles with eliminated x] on the (P, p) plane, where the graph of 
the equilibrium equation of state P(p) is also ploted. The unperturbed state corresponds to a point at the P(p) curve left to 
the upper stable point. Curves for t = 5, 10, 15, 20, 30 are shown. The short piece (10) of line near this point represents 
(splitting) P — p loci for times t = 5 and 10. 



FIG. 5. Same as in Fig. 3 for A = 30Aa, tr — 0.0033, uo = 0.03. Panel (a) shows velocity profiles for times t = 0, 2, 3, 4. 
Panel (b) shows density profiles for times t = 2, 3, 4. 



FIG. 6. P — p loci for the case shown in Fig. 5: curves are shown for t = 2, 2.5, 3, 4, 4.5, 5. 

FIG. 7. P c (t) for several values of f R ; plots show AP(t) = P c (t) - 1: f r = 0.33 (a), tr = 0.1 (b) and f R = 0.033 (c). At 
t < td Pc{t) monotonously decreases along P(p) curve to minimal value, at t ~ td it separates P(p) value and grows sharply to 
relax then to the isobaric distribution. Time is given in units ta- 
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